home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATHEMAT / STATISTI / 0910.ZIP / STAT-SAK.FOR < prev   
Text File  |  1987-02-11  |  56KB  |  1,902 lines
  1. C                                 STAT-SAK
  2. C                            The Statistician's
  3. C                             Swiss Army Knife   
  4. C                               Version 2.15
  5. C                            February 11, 1987
  6. C
  7. C                       "One of many STATOOLS(tm)..."
  8. C                                    by
  9. C
  10. C                             Gerard E. Dallal
  11. C                             53 Beltran Street
  12. C                             Malden, MA  02148
  13. C
  14. C
  15. C                                  NOTICE
  16. C
  17. C       Documentation and original code  copyright  1985  and 1986 by  
  18. C       Gerard E. Dallal. Reproduction of material for non-commercial
  19. C       purposes is permitted, without charge, provided that suitable 
  20. C       reference is made to STAT-SAK and its author.  
  21. C
  22. C       Neither  STAT-SAK nor its documentation should be modified in 
  23. C       any way without permission from the author,  except for those 
  24. C       changes  that  are  essential  to move   STAT-SAK to  another 
  25. C       computer.  
  26. C
  27.       CHARACTER*2 QUERY
  28.       DATA IIN /0/, IOUT /0/, NOPT /10/, QUERY /'  '/
  29. C
  30.       WRITE(IOUT,1)
  31.     1 FORMAT (////35X,'STAT-SAK'/
  32.      *   30X,'The Statistician''s'/
  33.      *   31X,'Swiss Army Knife'/
  34.      *   33X,'Version 2.15'/29X,'(c) 1985, 1986, 1987'//
  35.      *   26X,'"One of many STATOOLS(tm)..."'/
  36.      *   38X,'by'//
  37.      *   31X,'Gerard E. Dallal'/
  38.      *   31X,'53 Beltran Street'/
  39.      *   31X,'Malden, MA  02148'/////)
  40. C
  41.    30 WRITE(IOUT,40) 
  42.    40 FORMAT (/' Enter ''Q'' to quit, press  <ENTER>  to continue:  '$)
  43.       READ (IIN,50) QUERY
  44.    50 FORMAT (A2)
  45.    55 CALL UPCASE(QUERY)
  46.       IF (QUERY.EQ.'Q ') STOP ' '
  47.       IF (QUERY.EQ.'  ') GOTO 60
  48.       CALL OPTION(QUERY,NPICK)
  49.       
  50.       IF (NPICK.LT.1 .OR. NPICK.GT.NOPT) GOTO 60
  51. C
  52.       GOTO (100,200,300,400,500,600,700,800,900,1000), NPICK
  53. C
  54.  
  55. C
  56.    60 WRITE(IOUT,70)
  57.    70 FORMAT(/' Your options are:'//
  58.      *        '  DF -- distribution functions'/
  59.      *        '  CT -- independence in ',
  60.      *        'two-dimensional contingency tables'/
  61.      *        '  FI -- Fisher''s exact test for 2*2 tables'/
  62.      *        '  MH -- Mantel-Haenszel test / odds ratios'/
  63.      *        '  MC -- McNemar''s test'/
  64.      *        '  CC -- correlation coefficients'/
  65.      *        '  BP -- Bartholomew''s test for increasing proportions'/
  66.      *        '  BM -- Bartholomew''s test for increasing means'/
  67.      *        '  T1 -- one sample t-test from summary statistics'/
  68.      *        '  T2 -- two sample t-test from summary statistics'/)
  69. C
  70.       WRITE(IOUT,80)
  71.    80 FORMAT(' Your choice?:  ',$)
  72.       READ(IIN,50) QUERY
  73.       GOTO 55
  74. C
  75.   100 CALL PROB(IIN,IOUT)
  76.       GOTO 30
  77. C
  78.   200 CALL GOF(IIN,IOUT)
  79.       GOTO 30
  80. C
  81.   300 CALL MANTEL(IIN,IOUT)
  82.       GOTO 30
  83. C
  84.   400 CALL MCNEMR(IIN,IOUT)
  85.       GOTO 30 
  86. C
  87.   500 CALL CORR(IIN,IOUT)
  88.       GOTO 30
  89. C
  90.   600 CALL BART(IIN,IOUT)
  91.       GOTO 30
  92. C 
  93.   700 CALL TTEST1(IIN,IOUT)
  94.       GOTO 30
  95. C
  96.   800 CALL TTEST2(IIN,IOUT)
  97.       GOTO 30
  98. C
  99.   900 CALL OMEANS(IIN,IOUT)
  100.       GOTO 30
  101. C
  102.  1000 CALL FISHER(IIN,IOUT)
  103.       GOTO 30
  104. C
  105.       END
  106. C
  107.       SUBROUTINE UPCASE(TEXT)
  108.       CHARACTER*2 TEXT
  109.       CHARACTER*26 UPPER, LOWER
  110.       DATA UPPER / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' /
  111.       DATA LOWER / 'abcdefghijklmnopqrstuvwxyz' /
  112. C
  113.       DO 100 I = 1, 2
  114.       IF (TEXT(I:I).LE.'Z' .OR. TEXT(I:I).EQ.' ') GOTO 100
  115.       DO 90 K = 1, 26
  116.       IF (TEXT(I:I).NE.LOWER(K:K)) GOTO 90
  117.       TEXT(I:I) = UPPER(K:K)
  118.       GOTO 100
  119.    90 CONTINUE
  120.   100 CONTINUE
  121.       RETURN
  122.       END
  123. C
  124.       SUBROUTINE OPTION(TEXT,NO)
  125. C
  126.       CHARACTER*2 TEXT, LIST(10)
  127.       DATA NOPT /10/
  128. C      
  129.       DATA LIST(1), LIST(2), LIST(3), LIST(4)
  130.      *   /    'DF',    'CT',    'MH',    'MC'/
  131.       DATA LIST(5), LIST(6), LIST(7), LIST(8), LIST(9), LIST(10)
  132.      *  /     'CC',    'BP',    'T1',    'T2',    'BM',     'FI'/
  133. C
  134.       DO 100 I = 1, NOPT
  135.       NO = I
  136.       IF (TEXT.EQ.LIST(I)) RETURN
  137.   100 CONTINUE
  138.       NO = NO + 1
  139.       RETURN
  140.       END
  141. C
  142.       FUNCTION ALGAMA(S)
  143. C
  144. C        CACM ALGORITHM 291
  145. C        LOGARITHM OF GAMMA FUNCTION
  146. C        BY PIKE, M.C. AND HILL, I.D.
  147. C
  148. C        TRANSLATED INTO FORTRAN BY
  149. C                   Gerard E. Dallal
  150. C                   USDA-Human Nutrition Research Center on Aging
  151. C                              at Tufts University
  152. C                   711 Washington Street
  153. C                   Boston, MA  02111
  154. C
  155. C
  156.       X = S
  157.       ALGAMA = 0.0
  158.       IF (X .LE. 0.0) RETURN
  159.       F = 0.0
  160.       IF (X .GE. 7.0) GOTO 30
  161. C
  162.       F = 1.0
  163.       Z = X - 1.0
  164. C
  165.    10 Z = Z + 1.0
  166.       IF (Z .GE. 7.0) GOTO 20
  167.       X = Z
  168.       F = F * Z
  169.       GOTO 10
  170. C
  171.    20 X = X + 1
  172.       F = -ALOG(F)
  173. C
  174.    30 Z = 1.0 / X ** 2
  175.       ALGAMA = F + ( X - 0.5) * ALOG(X) - X + 0.918938533204673 +
  176.      1   (((-0.000595238095238 * Z + 0.000793650793651) * Z
  177.      2   - 0.002777777777778) * Z + 0.083333333333333) / X
  178.       RETURN
  179.       END
  180. C
  181.       FUNCTION ALNORM(X,UPPER)
  182. C
  183. C        ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, NO.3
  184. C
  185. C        EVALUATES THE TAIL AREA OF THE STANDARD NORMAL CURVE
  186. C        FROM X TO INFINITY IF UPPER IS .TRUE. OR
  187. C        FROM MINUS INFINITY TO X IF UPPER IS .FALSE.
  188. C
  189.       REAL LTONE, UTZERO, ZERO, HALF, ONE, CON, Z, Y, X
  190.       LOGICAL UPPER, UP
  191. C
  192. C        LTONE AND UTZERO MUST BE SET TO SUIT THE PARTICULAR
  193. C        COMPUTER (SEE INTRODUCTORY TEXT)
  194. C
  195.       DATA LTONE, UTZERO/7.0, 18.66/
  196.       DATA ZERO, HALF, ONE, CON/0.0, 0.5, 1.0, 1.28/
  197.       UP = UPPER
  198.       Z = X
  199.       IF (Z .GE. ZERO) GOTO 10
  200.       UP = .NOT. UP
  201.       Z = -Z
  202.    10 IF (Z .LE. LTONE .OR. UP .AND. Z .LE. UTZERO) GOTO 20
  203.       ALNORM = ZERO
  204.       GOTO 40
  205.    20 Y = HALF * Z * Z
  206.       IF (Z .GT. CON) GOTO 30
  207. C
  208.       ALNORM = HALF - Z * (0.398942280444 - 0.399903438504 * Y /
  209.      1   (Y + 5.75885480458 - 29.8213557808 /
  210.      2   (Y + 2.62433121679 + 48.6959930692 /
  211.      3   (Y + 5.92885724438))))
  212.       GOTO 40
  213. C
  214.    30 ALNORM = .398942280385 * EXP(-Y) /
  215.      1   (Z - 3.8052E-8 + 1.00000615302 /
  216.      2   (Z + 3.98064794E-4 + 1.98615381364 /
  217.      3   (Z - 0.151679116635 + 5.29330324926 /
  218.      4   (Z + 4.8385912808 - 15.1508972451 /
  219.      5   (Z + 0.742380924027 + 30.789933034 / (Z + 3.99019417011))))))
  220. C
  221.    40 IF (.NOT. UP) ALNORM = ONE - ALNORM
  222.       RETURN
  223.       END
  224. C
  225.       SUBROUTINE BART(IIN,IOUT)
  226. C
  227. C              BARTHOLOMEW'S TEST FOR ORDERED SAMPLES (INCREASING)
  228. C
  229.       DIMENSION CT(10),X(2,10),P(10)
  230.       CHARACTER*1 QUERY
  231. C
  232.     3 WRITE(IOUT,4)
  233.     4 FORMAT(/' Enter the number of columns [.LE.10]:  ',$)
  234.       READ(IIN,*,ERR=3) NCOL
  235.       IF(NCOL.LT.2 .OR. NCOL.GT.10) GOTO 3
  236. C
  237.     8 WRITE (IOUT,9)
  238.     9 FORMAT(' Do you wish to enter two rows (R) or ',
  239.      *       'one row with column totals (T)?:  '$)
  240.       READ (IIN,10) QUERY
  241.    10 FORMAT(A1)
  242.       IF (QUERY.NE.'R' .AND. QUERY.NE.'r' .AND. 
  243.      *   QUERY.NE.'T' .AND. QUERY.NE.'t') GO TO 8
  244. C
  245.       WRITE (IOUT,5)
  246.     5 FORMAT(' ')
  247.     6 WRITE(IOUT,7) 
  248.     7 FORMAT('         Enter row 1:  ',$)
  249.       READ(IIN,*,ERR=6) (X(1,J), J = 1, NCOL)
  250.       IF (QUERY.EQ.'2') GOTO 16
  251. C
  252.       IF (QUERY.EQ.'R' .OR. QUERY.EQ.'r') GOTO 16
  253.    11 WRITE (IOUT,12)
  254.    12 FORMAT(' Enter column totals:  ',$)
  255.       READ(IIN,*,ERR=11) (CT(J), J = 1, NCOL)
  256.       GOTO 18
  257. C
  258.    16 WRITE(IOUT,17) 
  259.    17 FORMAT('         Enter row 2:  ',$)
  260.       READ(IIN,*,ERR=16) (X(2,J), J = 1, NCOL)
  261. C
  262. C              GET COL TOTALS AND P'S
  263. C
  264.    18 PNUM=0.
  265.       PDEN=0.
  266.       DO 19 J=1,NCOL
  267.       IF (QUERY.EQ.'R' .OR. QUERY.EQ.'r') CT(J)=X(1,J)+X(2,J)
  268.       P(J)=X(1,J)/CT(J)
  269.       PNUM=PNUM+X(1,J)
  270.       PDEN=PDEN+CT(J)
  271.    19 CONTINUE
  272.       PBAR=PNUM/PDEN
  273. C
  274. C              PEARSON CHI-SQUARE STATISTIC
  275. C
  276.       PX2 = 0.0
  277.       DO 25 J = 1, NCOL
  278.       E11 = PBAR * CT(J)
  279.       PX2 = Px2 + (X(1,J) - E11)**2 * (1.0 / E11 + 1.0 / (CT(J) - E11))
  280.    25 CONTINUE
  281. C
  282. C              BARTHOLOMEW'S TEST
  283. C              FIRST MASSAGE TABLE
  284. C
  285.       DO 50 J=2,NCOL
  286.       IF (P(J-1).LE.P(J)) GO TO 50
  287. C
  288. C              P DECREASES...DETERMINE VALUE BY 'COLLAPSING'
  289. C              SO THAT SEQUENCE IS MONOTONE
  290. C
  291.       PNUM=P(J-1)*CT(J-1)+P(J)*CT(J)
  292.       PDEN=CT(J-1)+CT(J)
  293.       PX=PNUM/PDEN
  294. C
  295. C              IS IT NECESSARY TO GO FURTHER BACK?
  296. C
  297.       JS=J-1
  298.       IF (JS.EQ.1) GO TO 30
  299.       JM2=J-2
  300.       DO 29 I=1,JM2
  301.       IF (P(JS-1).LE.PX) GO TO 30
  302.       JS=JS-1
  303.       PNUM=PNUM+P(JS)*CT(JS)
  304.       PDEN=PDEN+CT(JS)
  305.       PX=PNUM/PDEN
  306.    29 CONTINUE
  307.    30 DO 40 J1=JS,J
  308.    40 P(J1)=PX
  309.    50 CONTINUE
  310. C
  311. C              CALCULATE STATISTIC
  312. C
  313.       CHISQ=0.
  314.       DO 70 J=1,NCOL
  315.    70 CHISQ=CHISQ+CT(J)*(P(J)-PBAR)**2
  316.       CHISQ=CHISQ/(PBAR*(1.-PBAR))
  317. C
  318.       WRITE (IOUT,80)CHISQ
  319.    80 FORMAT(/' Chi-square statistic for order restriction = ',G11.4)
  320. C
  321. C              INDICES FOR USE WITH TABLES
  322. C
  323.       CALL PBART(CHISQ,PVAL,NCOL,DIFF,2,CT)
  324.       IF (NCOL.EQ.3 .OR. NCOL.EQ.4) WRITE(IOUT,90) PVAL
  325.    90 FORMAT ('                            P-value (exact) = ',F7.3)
  326.       IF (NCOL.NE.3 .AND.NCOL.NE.4) WRITE(IOUT,100) PVAL
  327.   100 FORMAT ('           P-value (assuming equal weights) = ',F7.3)
  328. C
  329.   200 DIFF = PX2 - CHISQ
  330.       PDIFF = CHIDFC(DIFF,NCOL-1)
  331.       WRITE(IOUT,210) PX2, DIFF, PDIFF
  332.   210 FORMAT(/' Pearson chi-square statistic =',F8.2,5X,
  333.      * /'     Difference (lack of fit) =',F8.2,5X,'(P-value <',F6.3,')')
  334. C
  335.       RETURN
  336.       END
  337. C
  338.       FUNCTION BETAIN(X,P,Q)
  339. C
  340. C        ALGORITHM AS 63  APPL.STATIST. (1973), VOL.22, NO.3
  341. C        MODIFIED AS PER REMARK ASR 19  (1977), VOL.26, NO. 1
  342. C        ***[FAULT INDICATOR REMOVED BY G.E.D.]***
  343. C
  344. C        COMPUTES INCOMPLETE BETA FUNCTION RATIO FOR ARGUMENTS
  345. C        X BETWEEN ZERO AND ONE, P AND Q POSITIVE.
  346. C        LOG OF COMPLETE BETA FUNCTION, BETA, ASSUMED TO BE KNOWN.
  347. C        ***[CALCULATION OF BETA ADDED BY G.E.D]***
  348. C
  349.       LOGICAL INDEX
  350. C
  351. C        DEFINE ACCURACY AND INITIALIZE
  352. C
  353.       DATA ACU /0.1E-7/
  354.       BETAIN = X 
  355. C
  356. C        TEST FOR ADMISIBILITY OF AGRUMENTS
  357. C
  358.       IF (P .LE. 0.0 .OR. Q .LE. 0.0) STOP 40
  359.       IF (X .LT. 0.0 .OR .X .GT. 1.0) STOP 41
  360.       IF (X .EQ .0.0 .OR .X .EQ. 1.0) RETURN
  361. C
  362. C     CHANGE TAIL IF NECESSARY AND DETERMINE S
  363. C
  364.       PSQ = P + Q
  365.       CX = 1.0 - X
  366.       IF (P .GE .PSQ * X) GOTO 1
  367.       XX = CX
  368.       CX = X
  369.       PP = Q
  370.       QQ = P
  371.       INDEX = .TRUE.
  372.       GOTO 2
  373.     1 XX = X
  374.       PP = P
  375.       QQ = Q
  376.       INDEX = .FALSE.
  377.     2 TERM = 1.0
  378.       AI = 1.0
  379.       BETAIN = 1.0
  380.       NS = QQ + CX * PSQ
  381. C
  382. C        USE SOPER'S REDUCTION FORMULAE
  383. C
  384.       RX = XX / CX
  385.     3 TEMP = QQ - AI
  386.       IF (NS .EQ. 0) RX = XX
  387.     4 TERM = TERM * TEMP * RX / (PP + AI)
  388.       BETAIN = BETAIN + TERM
  389.       TEMP = ABS(TERM)
  390.       IF (TEMP .LE. ACU .AND. TEMP .LE. ACU * BETAIN) GOTO 5
  391.       AI = AI + 1.0
  392.       NS = NS - 1
  393.       IF (NS .GE. 0) GOTO 3
  394.       TEMP = PSQ
  395.       PSQ = PSQ + 1.0
  396.       GOTO 4
  397. C
  398. C        CALCULATE RESULT
  399. C
  400.     5 BETA = ALGAMA(P) + ALGAMA(Q) - ALGAMA(P+Q)
  401.       BETAIN = BETAIN * 
  402.      *     EXP(PP * ALOG(XX) + (QQ - 1.0) * ALOG(CX) - BETA) / PP
  403.       IF (INDEX) BETAIN = 1.0 - BETAIN
  404.       RETURN
  405.       END
  406. C
  407.       FUNCTION CHIDFC(X,NDF)
  408. C
  409. C        RETURNS P(CHISQ(NDF).GT.X), WHERE CHISQ(NDF) IS A CHI-SQUARE
  410. C        RANDOM VARIABLE WITH 'NDF' DEGREES OF FREEDOM.
  411. C
  412. C        IF  X.LE.0, RETURNS THE VALUE 1.0
  413. C
  414. C        UPPER TAIL RECURRENCE FORMULA...
  415. C           Q(X/NDF)=Q(X/NDF-2)+X**(NDF/2-1)*EXP(-X/2)/GAMMA(NDF/2),
  416. C              WHERE Q(X/NDF) IS THE UPPER TAIL, (X, INFINITY), OF THE
  417. C              CHI-SQUARE DISTRIBUTION WITH 'NDF' DEGREES OF FREEDOM.
  418. C        SEE HANDBOOK OF MATHEMATICAL FUNCTIONS,
  419. C        NBS,55(1964),26.4.8
  420. C
  421. C        USES WILSON-HILFERTY APPROXIMATION FOR NDF.GT.NDFBIG
  422. C        NBS HANDBOOK, 24.4.17
  423. C           (X/NDF)**(1./3.) IS APPROXIMATLEY NORMALLY DISTRIBUTED
  424. C           WITH MEAN 1.-2./(9.*NDF) AND VARIANCE 2./(9.*NDF)
  425. C
  426. C        REQUIRES FUNCTIONS ALNORM, ALGAMA
  427. C
  428. C        WRITTEN BY  Gerard E. Dallal
  429. C                    USDA-Human Nutrition Research Center on Aging
  430. C                               at Tufts University
  431. C                    711 Washington Street
  432. C                    Boston, MA  02111
  433. C
  434. C
  435. C        EXPBIG IS A NUMBER GREATER THAN THE MOST NEGATIVE VALID
  436. C           ARGUMENT TO THE EXP() FUNCTION.
  437. C
  438.       DATA NDFBIG /60/, EXPBIG /-80.0/
  439. C
  440.       CHIDFC = 1.0
  441.       IF (X .LE. 0.0 .OR. NDF .LT. 1) RETURN
  442.       CHIDFC = 0.
  443.       DF = NDF
  444.       IF (NDF .GT. NDFBIG) GOTO 50
  445.       HX = X / 2.0
  446.       IF ((NDF / 2) * 2 .EQ. NDF) GOTO 30
  447. C
  448. C        ODD NUMBER OF DEGREES OF FREEDOM
  449. C
  450.       IF (NDF .EQ. 1) GOTO 20
  451. C
  452.       INDEX = (NDF - 1) / 2
  453.       DO 10 I = 1, INDEX
  454.       C = FLOAT(I) + 0.5
  455.       D = (C - 1.0) * ALOG(HX) - HX - ALGAMA(C)
  456.       IF (D .GT. EXPBIG) CHIDFC = CHIDFC + EXP(D)
  457.    10 CONTINUE
  458.    20 CHIDFC = CHIDFC + 2.0 * ALNORM (SQRT(X), .TRUE.)
  459.       RETURN
  460. C
  461. C        EVEN NUMBER OF DEGREES OF FREEDOM
  462. C
  463.    30 INDEX = NDF / 2
  464.       DO 40 I = 1, INDEX
  465.       C = I
  466.       D =(C - 1.0) * ALOG(HX) - HX - ALGAMA(C)
  467.       IF (D .GT. EXPBIG) CHIDFC = CHIDFC + EXP(D)
  468.    40 CONTINUE
  469.       RETURN
  470. C
  471.    50 D = 2.0 / (9.0 * DF)
  472.       D = ((X / DF) ** (1.0 / 3.0) + D - 1.0) / SQRT(D)
  473.       CHIDFC = ALNORM(D, .TRUE.)
  474.       RETURN
  475.       END
  476. C
  477.       SUBROUTINE CORR(IIN,IOUT)
  478. C
  479.       CHARACTER*1 QUERY
  480. C
  481.    10 WRITE (IOUT,20)
  482.    20 FORMAT (/'   0 -- test that a population correlation '
  483.      *   'coefficient is zero'/
  484.      *         '   1 -- confidence interval for a single '
  485.      *   'correlation coefficient'/
  486.      *         '   2 -- compare two independent correlation '
  487.      *   'coefficients'//' Pick one:  ',$)
  488.       READ (IIN,30) QUERY
  489.    30 FORMAT (A1)
  490.       IF (QUERY.NE.'0' .AND. QUERY.NE.'1' .AND. QUERY.NE.'2') GOTO 10
  491. C
  492.       IF (QUERY.EQ.'0') CALL CORR0(IIN,IOUT)
  493.       IF (QUERY.EQ.'1') CALL CORR1(IIN,IOUT)
  494.       IF (QUERY.EQ.'2') CALL CORR2(IIN,IOUT)
  495.       RETURN
  496.       END
  497. C
  498.       SUBROUTINE CORR0(IIN,IOUT)
  499. C
  500.   110 WRITE(IOUT,120)
  501.   120 FORMAT(/' Enter sample correlation coefficient:  '$)
  502.       READ(IIN,*,ERR=110) R
  503.       IF(ABS(R).GT.1) GOTO 110
  504.   130 WRITE(IOUT,140)
  505.   140 FORMAT(' Enter number of observations:  '$)
  506.       READ(IIN,*,ERR=130) FN
  507.       IF (FN.LT.3.0) GOTO 130
  508. C
  509.       T = SQRT(FN - 2.0) * R / SQRT(1.0 - R**2)
  510.       NDF = FN - 2.0
  511.       P = FDFC(T**2,1.0,FLOAT(NDF))
  512.       WRITE (IOUT,210) T, NDF, P
  513.   210 FORMAT (
  514.      *   /'   t = ',G11.4,5X,'(',I4,' d.f.)',5X,'P-value = ',F6.4)
  515. C
  516.       RETURN
  517.       END
  518. C
  519.       SUBROUTINE CORR1(IIN,IOUT)
  520. C
  521. C        TWO-SIDED CONFIDENCE LIMITS FOR A POPULATION CORRELATION
  522. C        COEFFICIENT.  TRANSLATED FROM THE BASIC PROGRAM ON P.300 OF
  523. C          MAINDONALD, J.H.(1984)  STATISTICAL COMPUTATION.
  524. C          NEW YORK:  JOHN WILEY & SONS, INC.
  525. C
  526. C        USES APPROXIMATION FROM 
  527. C          WINTERBOTTOM, ALAN (1980).  ESTIMATION FOR THE BIVARIATE
  528. C          NORMAL CORRELATION COEFFICIENT USING ASYMPTOTIC EXPANSIONS 
  529. C
  530. C        THREE DIGIT ACCURACY FOR SAMPLES OF 10
  531. C        NEARLY FOUR DIGIT ACCURACY FOR SAMPLES OF 25 OR MORE     
  532. C
  533.       FNA(X,R,V,Z) = Z + X / SQRT(V) - R / (2.0 * V) + X * (X**2 +
  534.      *   3.0 * (1.0 + R**2)) / (12.0 * V * SQRT(V))
  535.       FNB(X,R,V) = R * (4.0 * (R * X)**2 + 5.0 * R**2 + 9.0) /
  536.      *   (24.0 * V**2)
  537.       FNC(X,R3,R4,V) = X * (X**4 + R3 * X**2 + R4) / 
  538.      *   (480.0 * V**2 * SQRT(V))
  539.       FNT(X,R,R3,R4,V,Z) = FNA(X,R,V,Z) - FNB(X,R,V) + FNC(X,R3,R4,V)
  540.       FNR(Z) = (EXP(2.0 * Z) - 1.0) / (EXP(2.0 * Z) + 1.0)
  541.  
  542. C
  543.    90 WRITE (IOUT,100)
  544.   100 FORMAT(/' Enter confidence level (0,1):  '$)
  545.       READ(IIN,*,ERR=90) CONF
  546.       IF (CONF.LE.0.0 .OR. CONF.GE.1.0) GOTO 90
  547.       Z0 = -GAUINV((1.0 - CONF)/ 2.0,IFAULT)
  548. C
  549.   110 WRITE(IOUT,120)
  550.   120 FORMAT(' Enter sample correlation coefficient'
  551.      *   'and number of observations:  '$)
  552.       READ(IIN,*,ERR=110) R,N2
  553.       IF(ABS(R).GE.1 .OR. N2.LE.2) GOTO 110
  554. C
  555.       V = N2 - 1
  556.       Z = 0.5 * ALOG((1.0 + R) / (1.0 - R))
  557.       R3 = 60.0 * R**4 - 30.0 * R**2 + 20.0
  558.       R4 = 165.0 * R**4 + 30.0 * R**2 + 15.0
  559.       X = -Z0
  560.       Z1 = FNT(X,R,R3,R4,V,Z)
  561.       X = Z0
  562.       Z2 = FNT(X,R,R3,R4,V,Z)
  563.       R1 = FNR(Z1)
  564.       R2 = FNR(Z2)
  565.       WRITE(IOUT,150) 100.0*CONF,R1,R2
  566.   150 FORMAT (/' The ',F4.1,'-% confidence interval is (',
  567.      *   F7.4,', ',F7.4,')')
  568.       RETURN
  569.       END
  570. C
  571.       SUBROUTINE CORR2(IIN,IOUT)
  572. C
  573. C             USES FISHER'S Z TO COMPARE TWO CORRELATION COEFFICENTS
  574. C
  575.    10 WRITE (IOUT,20) 
  576.    20 FORMAT(/' Enter first  correlation coefficient'
  577.      *   'and sample size:  '$)
  578.    30 READ (IIN,*,ERR=10) CC1, N1
  579.       IF(ABS(CC1).GE.1.0 .OR. N1.LE.2) GOTO 10
  580. C
  581.   110 WRITE (IOUT,120) 
  582.   120 FORMAT(/' Enter second correlation coefficient'
  583.      *   'and sample size:  '$)
  584.   130 READ (IIN,*,ERR=110) CC2, N2
  585.       IF(ABS(CC2).GE.1.0 .OR. N2.LE.2) GOTO 110
  586. C
  587. C
  588.       Z = 0.0
  589.       P = 0.0
  590.       IF (N1.LE.3 .OR. N2.LE.3) GOTO 400
  591. C
  592.       CC1 = 0.5 * ALOG((1.0 + CC1) / (1.0 - CC1))
  593.       CC2 = 0.5 * ALOG((1.0 + CC2) / (1.0 - CC2))
  594.       Z = ABS(CC1 - CC2)/ SQRT(1.0 / FLOAT(N1 - 3) + 
  595.      *   1.0 / FLOAT(N2 - 3))
  596.       P = 2.0 * ALNORM(Z,.TRUE.)
  597. C
  598.   400 WRITE(IOUT,410) Z, P
  599.   410 FORMAT (/' Test statistic = ',G11.4,'      P-value = ',F5.3)
  600.       RETURN
  601.       END
  602. C
  603.       FUNCTION FDFC(FQUAN, DFN, DFD)
  604. C
  605. C        THE COMPLEMENT OF THE CDF OF THE F DISTRIBUTION
  606. C        WITH DFN NUMERATOR DEGREES OF FREEDOM AND DFD DENOMINATOR 
  607. C        DEGREES OF FREEDOM EVALUATED AT FQUAN.
  608. C
  609. C
  610.       ESQ = DFD / (DFD + DFN * FQUAN)
  611.       FDFC = BETAIN (ESQ, DFD/2.0, DFN / 2.0)
  612.       RETURN
  613.       END
  614. C
  615.       SUBROUTINE FDFC0(IIN, IOUT, NPICK)
  616. C
  617. C        GET INFO FOR F-DISTRIBUTION CALCULATION
  618. C
  619.       IF(NPICK.EQ.1) GOTO 200
  620.    10 WRITE (IOUT,20)
  621.    20 FORMAT(/' Enter numerator degrees of freedom:  '$)
  622.       READ (IIN,*,ERR=10) DFN
  623.    30 WRITE (IOUT,40)
  624.    40 FORMAT(' Enter denominator degrees of freedom:  '$)
  625.       READ (IIN,*,ERR=30) DFD
  626.    50 WRITE (IOUT,60)
  627.    60 FORMAT(' Enter F quantile:  '$)
  628.       READ (IIN,*,ERR=50) FQUAN
  629.       PVAL = FDFC(FQUAN,DFN,DFD)
  630.       WRITE (IOUT,80) PVAL
  631.    80 FORMAT (/' Upper tail probability =',F8.5)
  632.       RETURN
  633. C
  634.   200 WRITE (IOUT,210)
  635.   210 FORMAT(/' Enter upper tail probability:  '$)
  636.       READ(IIN,*,ERR=200) PVAL
  637.       IF(PVAL.LT.0.0 .OR. PVAL.GT.1.0) GOTO 200
  638.   220 WRITE (IOUT,20)
  639.       READ (IIN,*,ERR=220) DFN
  640.   230 WRITE (IOUT,40)
  641.       READ (IIN,*,ERR=230) DFD
  642.       BETA = XINBTA(DFD/2.0,DFN/2.0,PVAL)
  643.       F = (1.0/BETA -1.0)*DFD/DFN
  644.       WRITE (IOUT,240) F
  645.   240 FORMAT(' Quantile:  ',G11.4)
  646.       RETURN   
  647.       END
  648. C
  649.       FUNCTION GAMAIN(X,P)
  650. C
  651. C        ALGORITHM AS 32 J.R.STATIST.SOC. C. (1970) VOL.19 NO.3
  652. C
  653. C        COMPUTES INCOMPLETE GAMMA RATIO FOR POSITIVE VALUES OF
  654. C        ARGUMENTS X AND P.  G MUST BE SUPPLIED AND SHOULD BE EQUAL TO
  655. C        LN(GAMMA(P)).
  656. C        IFAULT = 1 IF P.LE.0 ELSE 2 IF X.LT.0 ELSE 0
  657. C        USES SERIES EXPANSION IF P.GT.X OR X.LE.1, OTHERWISE A 
  658. C        CONTINUED FRACTION APPROXIMATION.C
  659. C
  660. C        [FAULT INDICATOR REMOVED BY G.E.D
  661. C        CALCULATION OF G INSERTED BY G.E.D.]
  662. C
  663.       DIMENSION PN(6)
  664. C
  665. C        DEFINE ACCURACY AND INITIALIZE
  666. C
  667.       ACU = 1.0E-8
  668.       OFLO = 1.0E30
  669.       GIN = 0.0
  670. C
  671. C        TEST FOR ADMISSIBILITY OF ARGUMENTS
  672. C
  673.       IF (P .LE. 0.0) STOP 21
  674.       IF (X .LT. 0.0) STOP 22
  675.       IF (X .EQ. 0.0) GOTO 50
  676.       G = ALGAMA(P)
  677.       FACTOR = EXP(P * ALOG(X) - X - G)
  678.       IF (X .GT. 1.0 .AND. X .GE. P) GOTO 30
  679. C
  680. C        CALCULATION BY SERIES EXPANSION
  681. C
  682.       GIN = 1.0
  683.       TERM = 1.0
  684.       RN = P
  685.    20 RN = RN + 1.0
  686.       TERM = TERM * X / RN
  687.       GIN = GIN + TERM
  688.       IF (TERM .GT. ACU) GOTO 20
  689.       GIN = GIN * FACTOR / P
  690.       GOTO 50
  691. C
  692. C        CALCULATION BY CONTINUED FRACTION
  693. C
  694.    30 A = 1.0 - P
  695.       B = A + X + 1.0
  696.       TERM = 0.0
  697.       PN(1) = 1.0
  698.       PN(2) = X
  699.       PN(3) = X + 1.0
  700.       PN(4) = X * B
  701.       GIN = PN(3) / PN(4)
  702.    32 A = A + 1.0
  703.       B = B + 2.0
  704.       TERM = TERM + 1.0
  705.       AN = A * TERM
  706.       DO 33 I = 1, 2
  707.    33 PN(I+4) = B * PN(I+2) - AN * PN(I)
  708.       IF (PN(6) .EQ. 0.0) GOTO 35
  709.       RN = PN(5) / PN(6)
  710.       DIF = ABS(GIN - RN)
  711.       IF (DIF .GT. ACU) GOTO 34
  712.       IF (DIF .LE. ACU * RN) GOTO 42
  713.    34 GIN = RN
  714.    35 DO 36 I = 1, 4
  715.    36 PN(I) = PN(I+2)
  716.       IF (ABS(PN(5)).LT.OFLO) GOTO 32
  717.       DO 41 I = 1, 4
  718.    41 PN(I) = PN(I) / OFLO
  719.       GOTO 32
  720.    42 GIN = 1.0 - FACTOR * GIN
  721. C
  722.    50 GAMAIN = GIN
  723.       RETURN
  724.       END
  725. C
  726.       FUNCTION GAUINV(P,IFAULT)
  727. C
  728. C              ALGORITHM AS 70 APPL. STATIST. (1974) VOL. 23, NO.1
  729. C              GAUINV FINDS PERCENTAGE POINTS OF THE NORMAL DISTRIBUTION
  730. C
  731.       DATA ZERO,ONE,HALF,ALIMIT/0.0, 1.0, 0.5, 1.0E-20/
  732. C
  733.       DATA             P0,  P1,             P2,                P3
  734.      1   / -.322232431088, -1., -.342242088547, -.204231210245E-1/
  735. C
  736.       DATA                P4,               Q0,            Q1
  737.      1   / -.453642210148E-4, .993484626060E-1, .588581570495/
  738. C
  739.       DATA            Q2,            Q3,              Q4
  740.      1   / .531103462366, .103537752850, .38560700634E-2/
  741. C
  742.       IFAULT=1
  743.       GAUINV=ZERO
  744.       PS=P
  745.       IF (PS.GT.HALF) PS=ONE-PS
  746.       IF (PS.LT.ALIMIT) RETURN
  747.       IFAULT=0
  748.       IF (PS.EQ.HALF) RETURN
  749.       YI=SQRT(ALOG(ONE/(PS*PS)))
  750.       GAUINV=YI+((((YI*P4+P3)*YI+P2)*YI+P1)*YI+P0)
  751.      1         /((((YI*Q4+Q3)*YI+Q2)*YI+Q1)*YI+Q0)
  752.       IF (P.LT.HALF) GAUINV=-GAUINV
  753.       RETURN
  754.       END
  755. C
  756.       SUBROUTINE GOF(IIN,IOUT)
  757. C
  758. C        TEST OF INDEPENDENCE AND HOMOGENEITY OF VARIANCE
  759. C        FOR TWO-DIMENSIONAL CONTINGENCY TABLES
  760. C
  761. C        LIMITATIONS:  NO MORE THAN FIVE ROWS OR COLUMNS.
  762. C
  763.       INTEGER ITABLE(5,10), IR(5), IC(10)
  764.       CHARACTER*1 QUERY
  765.       DATA NROW /5/, NCOL /10/
  766. C
  767.    30 WRITE(IOUT,70) NROW
  768.    70 FORMAT(/2X,'Enter the number of rows [.LE.',I1,']:  ',$)
  769.       READ(IIN,*,ERR=30) LROW
  770.    75 WRITE(IOUT,80) NCOL
  771.    80 FORMAT(2X,'Enter the number of columns [.LE.',I2,']:  ',$)
  772.       READ(IIN,*,ERR=75)LCOL
  773.       IF(LROW.LE.NROW.AND.LCOL.LE.NCOL)GO TO 100
  774.       WRITE(IOUT,90)
  775.    90 FORMAT(2X,'*** Permissible table dimensions exceeded')
  776.       GO TO 30
  777. C
  778.   100 WRITE (IOUT,110)
  779.   110 FORMAT(' ')
  780.       DO 130 I = 1, LROW
  781.   125 WRITE(IOUT,120) I
  782.   120 FORMAT(2X,'Enter row #',I1,':  ',$)
  783.       READ(IIN,*,ERR=125)(ITABLE(I,J), J = 1, LCOL)
  784.   130 CONTINUE
  785.  
  786.       DO 140 J = 1, LCOL
  787.       IC(J) = 0
  788.       DO 140 I = 1, LROW
  789.       IC(J) = IC(J) + ITABLE(I,J)
  790.   140 CONTINUE
  791. C
  792. C        COMPUTE ROW SUMS AND GRAND TOTAL
  793. C
  794.       GRAND = 0.0
  795.       DO 160 I = 1, LROW
  796.       IR(I) = 0
  797.       DO 150 J = 1, LCOL
  798.       IR(I) = IR(I) + ITABLE(I,J)
  799.   150 CONTINUE
  800.       GRAND = GRAND + FLOAT(IR(I))
  801.   160 CONTINUE
  802. C
  803. C        COMPUTE CHI-SQUARE P-VALUES
  804. C
  805.   180 CHISQ = 0.0
  806.       DO 190 I = 1, LROW
  807.       DO 190 J = 1, LCOL
  808.       EXPECT = FLOAT(IR(I)) * FLOAT(IC(J)) / GRAND
  809.       CHISQ = CHISQ + (FLOAT(ITABLE(I,J)) - EXPECT)**2 / EXPECT
  810.   190 CONTINUE
  811.       NDF = (LROW - 1) * (LCOL - 1)
  812.       PCHI = CHIDFC(CHISQ,NDF)
  813.       WRITE(IOUT,200) CHISQ,NDF,PCHI
  814.   200 FORMAT(/'  Pearson chi-square equals ',G11.4,' with',I3,
  815.      *   ' d.f.  (P-value = ',F7.5,')')
  816. C
  817.       IF (LROW.GT.2 .OR. LCOL.GT.2) GOTO 290
  818. C
  819. C        YATES'S CONTINUITY CORRECTED CHI-SQUARE STATISTIC
  820. C
  821.       YCHISQ = 0.0
  822.       DO 201 I = 1, 2
  823.       DO 201 J = 1, 2
  824.       EXPECT = FLOAT(IR(I)) * FLOAT(IC(J)) / GRAND
  825.       YCHISQ = YCHISQ + (ABS(FLOAT(ITABLE(I,J)) - EXPECT) - 0.5)**2 
  826.      *   / EXPECT
  827.   201 CONTINUE
  828.       YPCHI = CHIDFC(YCHISQ,1)
  829.       WRITE(IOUT,202) YCHISQ,YPCHI
  830.   202 FORMAT ('  Yates''s corrected chi-square equals ',G11.4,
  831.      *   '     (P-value = ',F7.5,')')
  832.  
  833.   290 WRITE(IOUT,300)
  834.   300 FORMAT(/2X,'Entering another table?  (Y or N):  '$)
  835.       READ(IIN,310) QUERY
  836.   310 FORMAT (A1)
  837.       IF (QUERY.EQ.'Y' .OR. QUERY.EQ.'y') GOTO 30
  838.       RETURN
  839.       END
  840. C
  841.       SUBROUTINE MCNEMR(IIN,IOUT)
  842. C
  843. C        MCNEMAR'S TEST 
  844. C        EXACT--USES BINOMIAL(n,0.5) DISTRIBUTION
  845. C
  846.    90 WRITE(IOUT,100)
  847.   100 FORMAT (/' Enter two discordant cells:  '$)
  848.       READ(IIN,*,ERR=90) B, C
  849. C
  850.       XNOBS = B + C
  851.       XNOBS1 = XNOBS + 1.0
  852.       NSTAT = MIN1(B,C)
  853.       PVAL = 0.0
  854.       DO 70 I = 0, NSTAT
  855.       TERM = ALGAMA(XNOBS1) - ALGAMA(FLOAT(I+1)) - 
  856.      *   ALGAMA(XNOBS1-FLOAT(I)) - XNOBS * ALOG(2.0)
  857.       IF (TERM.GT.-78.0) PVAL = PVAL + EXP(TERM)
  858.    70 CONTINUE
  859.       PVAL = AMIN1(1.0, 2.0 * PVAL)
  860. C
  861.       WRITE(IOUT,120) PVAL
  862.   120 FORMAT(/' P-value = ',F6.4)
  863.       RETURN
  864.       END
  865. C
  866.       FUNCTION PPCHI2(P,V)
  867. C
  868. C        ALGORITHM AS 91  APPL. STATIST. (1975) VOL.24, NO.3
  869. C
  870. C        TO EVALUATE THE PECENTAGE POINTS OF THE CHI-SQUARED
  871. C        PROBABILITY DISTRIBUTION FUNCTION.
  872. C        P MUST LIE IN THE RANGE 0.000002 TO 0.999998, V MUST BE POSITIVE,
  873. C        G MUST BE SUPPLIED AND SHOULD BE EQUAL TO LN(GAMMA(V/2.0))
  874. C
  875. C        [FAULT INDICATOR REMOVED BY G.E.D.  
  876. C        CALCULATION OF G ADDED BY G.E.D.
  877. C
  878.       DATA E, AA /0.5E-6, 0.6931471805/
  879. C
  880. C        AFTER DEFINING THE ACCURACY AND LN(2), TEST ARGUMENTS AND INITIALIZE
  881. C
  882.       PPCHI2 = -1.0
  883.       IF (P .LT. 0.000002 .OR. P .GT. 0.999998) STOP 31
  884.       IF (V .LE. 0.0) STOP 32
  885.       G = ALGAMA(V / 2.0)
  886.       XX = 0.5 * V
  887.       C = XX - 1.0
  888. C
  889. C        STARTING APPROXIMATION FOR SMALL CHI-SQUARED
  890. C
  891.       IF (V .GE. -1.24 * ALOG(P)) GOTO 1
  892.       CH = (P * XX * EXP(G + XX * AA)) ** (1.0 / XX)
  893.       IF (CH - E) 6, 4, 4
  894. C
  895. C        STARTING APPROXIMATION FOR V LESS THAN OR EQUAL TO 0.32
  896. C
  897.     1 IF (V .GT. 0.32) GOTO 3
  898.       CH = 0.4
  899.       A = ALOG(1.0 - P)
  900.     2 Q = CH
  901.       P1 = 1.0 + CH * (4.67 + CH)
  902.       P2 = CH * (6.73 + CH * (6.66 + CH))
  903.       T = -0.5 + (4.67 + 2.0 * CH) / P1 -
  904.      *   (6.73 + CH * (13.32 + 3.0 * CH)) / P2
  905.       CH = CH - (1.0 - EXP(A + G + 0.5 * CH + C * AA) * P2 / P1) / T
  906.       IF (ABS(Q / CH - 1.0) - 0.01) 4, 4, 2
  907. C
  908. C        CALL TO ALGORITHM AS 70 - NOTE THAT P HAS BEEN TESTED ABOVE
  909. C
  910.     3 X = GAUINV(P,IF1)
  911. C
  912. C        STARTING APPROXIMATION USING WILSON AND HILFERTY ESTIMATE
  913. C
  914.       P1 = 0.222222 / V
  915.       CH = V * (X * SQRT(P1) + 1.0 - P1) ** 3
  916. C
  917. C        STARTING APPROXIMATION FOR P TENDING TO 1
  918. C
  919.       IF (CH .GT. 2.2 * V + 6.0)
  920.      *   CH = -2.0 * (ALOG(1.0 - P) - C * ALOG(0.5 * CH) + G)
  921. C
  922. C        CALL TO ALGORITHM AS 32 AND CALCULATION OF SEVEN TERM
  923. C        TAYLOR SERIES
  924. C
  925.     4 Q = CH
  926.       P1 = 0.5 * CH
  927.       P2 = P - GAMAIN(P1, XX)
  928.       T = P2 * EXP(XX * AA + G + P1 - C * ALOG(CH))
  929.       B = T / CH
  930.       A = 0.5 * T - B * C
  931.       S1 = (210.0+A*(140.0+A*(105.0+A*(84.0+A*(70.0+60.0*A))))) / 420.0
  932.       S2 = (420.0+A*(735.0+A*(966.0+A*(1141.0+1278.0*A)))) / 2520.0
  933.       S3 = (210.0 + A * (426.0 + A * (707.0 + 932.0 * A))) / 2520.0
  934.       S4 =(252.0+A*(672.0+1182.0*A)+C*(294.0+A*(889.0+1740.0*A)))/5040.0
  935.       S5 = (84.0 + 264.0 * A + C * (175.0 + 606.0 * A)) / 2520.0
  936.       S6 = (120.0 + C * (346.0 + 127.0 * C)) / 5040.0
  937.       CH = CH+T*(1.0+0.5*T*S1-B*C*(S1-B*(S2-B*(S3-B*(S4-B*(S5-B*S6))))))
  938.       IF (ABS(Q / CH - 1.0) .GT. E) GOTO 4
  939. C
  940.     6 PPCHI2 = CH
  941.       RETURN
  942.       END
  943. C
  944.       SUBROUTINE PROB(IIN, IOUT)
  945. C
  946. C        PROBABILITY CALCULATION
  947. C
  948.       CHARACTER*2 QUERY
  949.       DATA MAXDIST /10/, KNORM /1/, KT /3/, KCHISQ /5/, KF /7/
  950. C
  951.       QUERY = ' '
  952.    10 WRITE(IOUT,20) 
  953.    20 FORMAT (/' Enter ''R'' to return to main menu,',
  954.      *   '  press  <ENTER>  to continue:  '$)
  955.    30 READ (IIN,40) QUERY
  956.    40 FORMAT (A2)
  957.       CALL UPCASE(QUERY)
  958.       IF (QUERY.EQ.'R ') RETURN
  959.       IF (QUERY.EQ.'Q ') STOP ' '
  960.       IF (QUERY.EQ.'  ') GOTO 50
  961.       CALL DFOPT(QUERY,NDIST)
  962.       IF (NPICK.GT.MAXDIST) GOTO 10
  963.       GOTO (100,100,200,200,300,300,400,400,500,600), NDIST
  964. C
  965.    50 WRITE (IOUT,60)
  966.    60 FORMAT (/' Your options are:'//
  967.      *   '       NQ -- standard normal:  quantile to probability'/
  968.      *   '       NP -- standard normal:  probability to quantile'/
  969.      *   '       TQ -- Student''s t:  quantile to probability'/
  970.      *   '       TP -- Student''s t:  probability to quantile'/
  971.      *   '       CQ -- chi-square:  quantile to probability'/
  972.      *   '       CP -- chi-square:  probability to quantile'/
  973.      *   '       FQ -- F:  quantile to probability'/
  974.      *   '       FP -- F:  probability to quantile'/
  975.      *   '       BI -- binomial'/
  976.      *   '       P  -- Poisson'//
  977.      *   ' Your choice?:  '$)
  978.       GOTO 30
  979. C
  980.   100 CALL N01(IIN,IOUT,NDIST-KNORM)
  981.       GOTO 10
  982. C
  983.   200 CALL STUDNT(IIN,IOUT,NDIST-KT)
  984.       GOTO 10
  985. C
  986.   300 CALL X2DFC0(IIN,IOUT,NDIST-KCHISQ)
  987.       GOTO 10
  988. C
  989.   400 CALL FDFC0(IIN,IOUT,NDIST-KF)
  990.       GOTO 10
  991. C
  992.   500 CALL BINOM(IIN,IOUT)
  993.       GOTO 10
  994. C
  995.   600 CALL POISN(IIN,IOUT)
  996.       GOTO 10
  997.       END
  998. C
  999.       SUBROUTINE DFOPT(TEXT,NO)
  1000. C
  1001.       CHARACTER*2 TEXT, LIST(10)
  1002.       DATA NOPT /10/
  1003. C      
  1004.       DATA LIST(1), LIST(2), LIST(3), LIST(4), LIST(5)
  1005.      *   /    'NQ',    'NP',    'TQ',    'TP',    'CQ'/
  1006.       DATA LIST(6), LIST(7), LIST(8), LIST(9), LIST(10)
  1007.      *  /     'CP',    'FQ ',   'FP',    'BI',     'P '/
  1008. C
  1009.       DO 100 I = 1, NOPT
  1010.       NO = I
  1011.       IF (TEXT.EQ.LIST(I)) RETURN
  1012.   100 CONTINUE
  1013.       NO = NO + 1
  1014.       RETURN
  1015.       END
  1016. C
  1017.       SUBROUTINE TTEST1(IIN,IOUT)
  1018. C
  1019. C        ONE SAMPLE T TEST FROM SUMMARY STATISTICS
  1020. C
  1021.       CHARACTER*1 QUERY
  1022. C
  1023.    10 WRITE (IOUT,20)
  1024.    20 FORMAT(/' Does the summary include a standard deviation (D) or '/
  1025.      *        '                              a standard error (E)?:  '$)
  1026.       READ (IIN,30) QUERY
  1027.    30 FORMAT(A1)
  1028.       IF (QUERY.NE.'D' .AND. QUERY.NE.'d' .AND. 
  1029.      *   QUERY.NE.'E' .AND. QUERY.NE.'e') GO TO 10
  1030. C
  1031.       IF (QUERY.EQ.'E' .OR. QUERY.EQ.'e') GOTO 140
  1032. C
  1033. C        STANDARD DEVIATION
  1034. C
  1035.    40 WRITE (IOUT,50)
  1036.    50 FORMAT (/' Enter the mean, standard deviation, and sample size:')
  1037.       READ (IIN,*,ERR=40) XBAR, SD, FN
  1038.       SE = SD / SQRT(FN)
  1039.       GOTO 200
  1040. C
  1041. C        STANDARD ERROR
  1042. C
  1043.   140 WRITE (IOUT,150)
  1044.   150 FORMAT (/' Enter the mean, standard error, and sample size:')
  1045.       READ (IIN,*,ERR=140) XBAR, SE, FN
  1046. C
  1047.   200 T = XBAR / SE
  1048. C
  1049.       NDF = FN - 1.0
  1050. C
  1051.       P = FDFC(T**2,1.0,FLOAT(NDF))
  1052. C
  1053.       WRITE (IOUT,210) T, NDF, P
  1054.   210 FORMAT (
  1055.      *   /'   t = ',G11.4,5X,'(',I4,' d.f.)',5X,'P-value = ',F6.4)
  1056.       RETURN
  1057.       END
  1058. C
  1059.       SUBROUTINE TTEST2(IIN,IOUT)
  1060. C
  1061. C        T TEST FOR INDEPENDENT SAMPLES FROM SUMMARY STATISTICS
  1062. C
  1063.       CHARACTER*1 QUERY
  1064. C
  1065.    10 WRITE (IOUT,20)
  1066.    20 FORMAT(/' Do the summaries include standard deviations (D) or '/
  1067.      *        '                              standard errors (E)?:  '$)
  1068.       READ (IIN,30) QUERY
  1069.    30 FORMAT(A1)
  1070.       IF (QUERY.NE.'D' .AND. QUERY.NE.'d' .AND. 
  1071.      *   QUERY.NE.'E' .AND. QUERY.NE.'e') GO TO 10
  1072. C
  1073.       IF (QUERY.EQ.'E' .OR. QUERY.EQ.'e') GOTO 140
  1074. C
  1075. C        STANDARD DEVIATIONS
  1076. C
  1077.    40 WRITE (IOUT,50)
  1078.    50 FORMAT (/' Enter the mean, standard deviation, and sample size',
  1079.      *   ' for sample 1:')
  1080.       READ (IIN,*,ERR=40) XBAR1, SD1, FN1
  1081.       SE1 = SD1 / SQRT(FN1)
  1082.    60 WRITE (IOUT,70)
  1083.    70 FORMAT (' Enter the mean, standard deviation, and sample size',
  1084.      *   ' for sample 2:')
  1085.       READ (IIN,*,ERR=60) XBAR2, SD2, FN2
  1086.       SE2 = SD2 / SQRT(FN2)
  1087.       GOTO 200
  1088. C
  1089. C        STANDARD ERRORS
  1090. C
  1091.   140 WRITE (IOUT,150)
  1092.   150 FORMAT (/' Enter the mean, standard error, and sample size',
  1093.      *   ' for sample 1:')
  1094.       READ (IIN,*,ERR=140) XBAR1, SE1, FN1
  1095.       SD1 = SE1 * SQRT(FN1)
  1096.   160 WRITE (IOUT,170)
  1097.   170 FORMAT (' Enter the mean, standard error, and sample size',
  1098.      *   ' for sample 2:')
  1099.       READ (IIN,*,ERR=160) XBAR2, SE2, FN2
  1100.       SD2 = SE2 * SQRT(FN2)
  1101. C
  1102. C        EQUAL VARIANCES
  1103. C
  1104.   200 XBARD = XBAR1 - XBAR2
  1105.       IFN1 = FN1
  1106.       IFN2 = FN2
  1107.       SDP = SQRT((((FN1 - 1.0) * SD1 * SD1 + (FN2 - 1.0) * SD2 * SD2)
  1108.      *   / (FN1 + FN2 - 2.0)))
  1109.       SEDEQ = SDP * SQRT(1.0 / FN1 + 1.0 / FN2)
  1110.       TEQ = XBARD / SEDEQ
  1111.       NDFEQ = FN1 + FN2 - 2.0
  1112.       PEQ = FDFC(TEQ**2,1.0,FLOAT(NDFEQ))
  1113. C
  1114.        WRITE(IOUT,210) XBAR1, XBAR2, XBARD, SD1, SD2, IFN1, IFN2
  1115.   210 FORMAT (/24X,'sample 1      sample 2     difference'/
  1116.      *   17X,'mean',3G14.5/
  1117.      *   3X,'standard deviation',2G14.5/
  1118.      *   10X,'sample size',I10,I14)
  1119. C
  1120. C        UNEQUAL VARIANCES
  1121. C
  1122.       IF (FN1.LE.1.0 .OR. FN2.LE.1.0 .OR. SD1.EQ.0.0 .OR. SD2.EQ.0.0) 
  1123.      *   GOTO 225
  1124.       SEDNEQ = SQRT(SE1 * SE1 + SE2 * SE2)
  1125.       TNEQ = XBARD / SEDNEQ
  1126.       C = SE1**2 / (SE1**2 + SE2**2)
  1127.       DFNEQ = 1.0 / ( C**2 / (FN1 - 1.0) + (1.0 - C)**2 / (FN2 - 1.0))
  1128.       PNEQ = FDFC(TNEQ**2, 1.0, DFNEQ)
  1129.       FRAT = SD1**2 / SD2**2
  1130.       NDFN = FN1 - 1.0
  1131.       NDFD = FN2 - 1.0
  1132.       IF (FRAT.GE.1.0) GOTO 219
  1133.       NDFD = FN1 - 1.0
  1134.       NDFN = FN2 - 1.0
  1135.       FRAT = 1.0 / FRAT
  1136.   219 PF =  2.0 * FDFC(FRAT,FLOAT(NDFN),FLOAT(NDFD))
  1137.       WRITE (IOUT,220) SDP, FRAT, NDFN, NDFD, PF
  1138.   220 FORMAT (/'      Pooled standard deviation:  ',G12.5//
  1139.      *   ' F-test for equality of variances:  F-ratio:  ',F7.2/
  1140.      *   '                             numerator d.f.:',I9/
  1141.      *   '                           denominator d.f.:',I9/
  1142.      *   '                                    P-value:  ',F7.4/)
  1143. C
  1144.   225 WRITE (IOUT,230) TEQ, NDFEQ, PEQ
  1145.   230 FORMAT (
  1146.      *   /'   Equal variances:  t = ',G11.4,5X,'(',I4,'   d.f.)',5X,
  1147.      *   'P-value = ',F6.4)
  1148.       IF (FN1.GT.1.0 .AND. FN2.GT.1.0 .AND. SD1.GT.0.0 .AND. 
  1149.      *   SD2.GT.0.0)  WRITE (IOUT,240) TNEQ, DFNEQ, PNEQ
  1150.   240 FORMAT (
  1151.      *    ' Unequal variances:  t = ',G11.4,5X,'(',F6.1,' d.f.)',5X,
  1152.      *   'P-value = ',F6.4)
  1153. C
  1154.       CIU = XINBTA(FLOAT(NDFEQ)/2.0, 0.5, 0.05)
  1155.       CIU = SQRT(FLOAT(NDFEQ)*(1.0/ CIU - 1.0))
  1156.       CIL = XBARD - CIU * SEDEQ
  1157.       CIU = XBARD + CIU * SEDEQ
  1158.       WRITE(IOUT,350) SEDEQ, CIL, CIU
  1159.   350 FORMAT (/'   Equal variances:  standard error of difference:  ',
  1160.      *   G12.5/
  1161.      *   21X,'95% C.I. for mean difference:  (',G12.5,' ,',G12.5')')
  1162.       IF (FN1.LE.1.0 .OR. FN2.LE.1.0 .OR. SD1.EQ.0.0 .OR. SD2.EQ.0.0) 
  1163.      *   RETURN
  1164. C
  1165.       CIU = XINBTA(DFNEQ/2.0, 0.5, 0.05)
  1166.       CIU = SQRT(DFNEQ*(1.0/ CIU - 1.0))
  1167.       CIL = XBARD - CIU * SEDNEQ
  1168.       CIU = XBARD + CIU * SEDNEQ
  1169.       WRITE(IOUT,360) SEDNEQ, CIL, CIU
  1170.   360 FORMAT (/' Unequal variances:  standard error of difference:  ',
  1171.      *   G12.5/
  1172.      *   21X,'95% C.I. for mean difference:  (',G12.5,' ,',G12.5')'/)
  1173. C
  1174.       RETURN
  1175.       END
  1176. C
  1177.       SUBROUTINE X2DFC0(IIN, IOUT, NPICK)
  1178. C
  1179. C        CHI-SQUARE DISTRIBUTION
  1180. C
  1181.       IF(NPICK.EQ.1) GOTO 200
  1182. C
  1183.    30 WRITE (IOUT,40)
  1184.    40 FORMAT(/' Enter chi-square quantile:  '$)
  1185.       READ (IIN,*,ERR=30) CHISQ
  1186. C
  1187.   110 WRITE (IOUT,120)
  1188.   120 FORMAT(' Enter degrees of freedom:  '$)
  1189.       READ (IIN,*,ERR=110) DF
  1190.       IF (DF.NE.FLOAT(IFIX(DF))) GOTO 110
  1191.       NDF = DF
  1192.       PVAL = CHIDFC(CHISQ,NDF)
  1193.       WRITE (IOUT,150) PVAL
  1194.   150 FORMAT (' Upper tail probability =',F8.5)
  1195.       RETURN
  1196. C
  1197.   200 WRITE (IOUT,210)
  1198.   210 FORMAT(/' Enter upper tail probability:  '$)
  1199.       READ(IIN,*,ERR=200) PVAL
  1200.       IF(PVAL.LT.0.0 .OR. PVAL.GT.1.0) GOTO 200
  1201.   220 WRITE (IOUT,230)
  1202.   230 FORMAT(' Enter degrees of freedom:  '$)
  1203.       READ (IIN,*,ERR=220) DF
  1204.       WRITE(IOUT,250) PPCHI2(1.0 - PVAL,DF)
  1205.   250 FORMAT (' Upper tail percentile = ',G12.5)
  1206.       RETURN
  1207.       END
  1208. C
  1209.       SUBROUTINE STUDNT(IIN, IOUT, NPICK)
  1210. C
  1211. C        INFORMATION FOR STUDENT'S T-DISTRIBUTION
  1212. C
  1213.       IF(NPICK.EQ.0) GOTO 230
  1214. C
  1215.    10 WRITE (IOUT,20)
  1216.    20 FORMAT(/' Enter upper tail probability:  '$)
  1217.       READ(IIN,*,ERR=10) PVAL
  1218.       IF(PVAL.LE.0.0 .OR. PVAL.GE.1.0) GOTO 10
  1219.    30 WRITE (IOUT,40)
  1220.    40 FORMAT(' Enter degrees of freedom:  '$)
  1221.       READ (IIN,*,ERR=30) DF
  1222.       T = XINBTA(DF/2.0, 0.5, 2.0 * AMIN1(PVAL,1.0-PVAL))
  1223.       T = SQRT(DF*(1.0/ T - 1.0))
  1224.       IF (PVAL.GT.0.5) T = -T
  1225.       WRITE(IOUT,50) T
  1226.    50 FORMAT (' Upper tail percentile: ',F8.3)
  1227.       RETURN
  1228. C
  1229.   230 WRITE (IOUT,240)
  1230.   240 FORMAT(/' Enter t quantile:  '$)
  1231.       READ (IIN,*,ERR=230) T
  1232.   250 WRITE (IOUT,260)
  1233.   260 FORMAT(' Enter degrees of freedom:  '$)
  1234.       READ (IIN,*,ERR=250) DF
  1235.       IF (DF.NE.FLOAT(IFIX(DF))) GOTO 250
  1236.       NDF = DF
  1237.       PVAL = 0.5 * FDFC(T**2,1.0,DF)
  1238.       IF (T.GT.0) PVAL = 1.0 - PVAL
  1239. C
  1240.       WRITE(IOUT,270) 1.0 - PVAL, PVAL, 
  1241.      *   2.0 * AMIN1(1.0 - PVAL, PVAL)
  1242.   270 FORMAT(/25X,'Upper tail = ',F7.4,
  1243.      *       /25X,'Lower tail = ',F7.4,
  1244.      *       /25X,'Two-tailed = ',F7.4)
  1245.       RETURN
  1246.       END
  1247. C
  1248.       SUBROUTINE N01(IIN, IOUT, NPICK)
  1249. C
  1250. C        INFORMATION FOR STANDARD NORMAL DISTRIBUTION
  1251. C
  1252.       IF(NPICK.EQ.1) GOTO 200
  1253. C
  1254.   130 WRITE (IOUT,140)
  1255.   140 FORMAT(/' Enter quantile:  '$)
  1256.       READ (IIN,*,ERR=130) Z
  1257.       ZC = Z      
  1258.       IF (Z.NE.0.0) GOTO 160
  1259.   150 WRITE (IOUT,155)
  1260.   155 FORMAT (' Enter A,B,C.  STAT-SAK will use A / (B / SQRT(C)):')
  1261.       READ (IIN,*,ERR=150) A,B,C
  1262.       Z = A / (B / SQRT(C))
  1263.   160 PVAL = ALNORM(Z,.TRUE.)
  1264. C
  1265.       IF (ZC.EQ.0.0) WRITE (IOUT,165) Z
  1266.   165 FORMAT(
  1267.      *       /25X,'         z = ',F7.3)
  1268.       WRITE(IOUT,170) PVAL, 1.0 - PVAL, 
  1269.      *   2.0 * AMIN1(1.0 - PVAL, PVAL)
  1270.   170 FORMAT(/25X,'Upper tail = ',F7.4,
  1271.      *       /25X,'Lower tail = ',F7.4,
  1272.      *       /25X,'Two-tailed = ',F7.4)
  1273.       RETURN
  1274. C
  1275.   200 WRITE (IOUT,210)
  1276.   210 FORMAT(/' Enter upper tail probability:  '$)
  1277.       READ(IIN,*,ERR=200) PVAL
  1278.       IF(PVAL.LT.0.0 .OR. PVAL.GT.1.0) GOTO 200
  1279.       WRITE(IOUT,250) ABS(GAUINV(PVAL,IFAULT))
  1280.   250 FORMAT (/' Quantile = ',G12.5)
  1281.       RETURN
  1282.       END
  1283. C
  1284.        SUBROUTINE OMEANS(IIN,IOUT)
  1285. C
  1286. C              BARTHOLOMEW'S TEST FOR ORDERED SAMPLES (INCREASING)
  1287. C
  1288.       DIMENSION XBAR(10), SD(10), SE(10), FN(10), DF(10)
  1289.       LOGICAL QSD
  1290.       CHARACTER*1 QUERY
  1291.       CHARACTER FNAME*50
  1292.       DATA IFIN /10/
  1293. C
  1294.     3 WRITE(IOUT,4)
  1295.     4 FORMAT(/' Enter the number of samples [.LE.10]:  ',$)
  1296.       READ(IIN,*,ERR=3) NSAM
  1297.       IF(NSAM.LT.2 .OR. NSAM.GT.10) GOTO 3
  1298.       TDF = NSAM - 1
  1299. C
  1300.     8 WRITE (IOUT,9)
  1301.     9 FORMAT(' Do you wish to enter standard errors (E) or ',
  1302.      *       'standard deviations (D)?:  '$)
  1303.       READ (IIN,10) QUERY
  1304.    10 FORMAT(A1)
  1305.       IF (QUERY.NE.'D' .AND. QUERY.NE.'d' .AND. 
  1306.      *   QUERY.NE.'E' .AND. QUERY.NE.'e') GO TO 8
  1307.       QSD = .TRUE.
  1308.       IF (QUERY.EQ.'E' .OR. QUERY.EQ.'e') QSD = .FALSE.
  1309. C
  1310.       WRITE (IOUT,900)
  1311.   900 FORMAT(' Are data contained in an external file?  (Y or N):  '$)
  1312.       READ (IIN,10) QUERY
  1313.       IF (QUERY.NE.'Y' .AND. QUERY.NE.'y') GOTO 999
  1314.   905 WRITE (IOUT,910)
  1315.   910 FORMAT (' Enter filename:  '$)
  1316.       READ (IIN,'(A)') FNAME
  1317.       OPEN (IFIN, FILE = FNAME, STATUS = 'OLD', IOSTAT =IOCHK)
  1318.       IF (IOCHK.NE.0) GOTO 905
  1319.       IF (QSD) GOTO 950
  1320. C
  1321.       WRITE (IOUT,920)
  1322.   920 FORMAT (/'    Mean      S.E.        count'/)
  1323.       DO 940 I = 1, NSAM
  1324.       READ (IFIN,*) XBAR(I), SE(I), FN(I)
  1325.       WRITE (IOUT,930) XBAR(I), SE(I), FN(I)
  1326.       SD(I) = SE(I) * SQRT(FN(I))
  1327.   930 FORMAT (1X,2G11.4,F8.0)
  1328.   940 CONTINUE
  1329.       CLOSE (IFIN)
  1330.       GOTO 26
  1331. C
  1332.   950 WRITE (IOUT,960)
  1333.   960 FORMAT (/'    Mean      S.D.        count'/)
  1334.       DO 970 I = 1, NSAM
  1335.       READ (IFIN,*) XBAR(I), SD(I), FN(I)
  1336.       WRITE (IOUT,930) XBAR(I), SD(I), FN(I)
  1337.       SE(I) = SD(I) / SQRT(FN(I))
  1338.   970 CONTINUE
  1339.       CLOSE (IFIN)
  1340.       GOTO 26      
  1341. C
  1342. C
  1343.   999 IF (QSD) GOTO 16
  1344.       DO 14 I = 1, NSAM
  1345.    11 WRITE (IOUT,12) I
  1346.    12 FORMAT(' Enter mean, standard error, and sample size for sample',
  1347.      *   I3,':')
  1348.       READ(IIN,*,ERR=11) XBAR(I), SE(I), FN(I)
  1349.       SD(I) = SE(I) * SQRT(FN(I))
  1350.    14 CONTINUE
  1351.       GOTO 26
  1352. C
  1353.    16 DO 25 I = 1, NSAM
  1354.    18 WRITE (IOUT,19) I
  1355.    19 FORMAT(' Enter mean, standard deviation, and sample size for ',
  1356.      *'sample',I3,':')
  1357.       READ(IIN,*,ERR=18) XBAR(I), SD(I), FN(I)
  1358.       SE(I) = SD(I) / SQRT(FN(I))
  1359.    25 CONTINUE
  1360. C
  1361. C              GET BETWEEN GROUP SUM OF SQUARES
  1362. C
  1363.    26 TBAR = 0.0
  1364.       TN = 0.0
  1365.       DO 27 I = 1, NSAM
  1366.       DF(I) = 1.0
  1367.       TBAR = TBAR + XBAR(I) * FN(I)
  1368.       TN = TN + FN(I)
  1369.    27 CONTINUE
  1370.       TBAR = TBAR / TN
  1371. C
  1372.       BSS = 0.0
  1373.       WSS = 0.0
  1374.       DO 28 I = 1, NSAM
  1375.       BSS = BSS + FN(I) * (XBAR(I) - TBAR)**2
  1376.       WSS = WSS + (FN(I) - 1.0) * SD(I)**2
  1377.    28 CONTINUE
  1378. C
  1379. C              BARTHOLOMEW'S TEST
  1380. C
  1381.       DO 50 J=2,NSAM
  1382.       IF (XBAR(J-1).LE.XBAR(J)) GO TO 50
  1383. C
  1384. C              MEAN DECREASES...DETERMINE VALUE BY POOLING
  1385. C              SO THAT SEQUENCE IS MONOTONE
  1386. C
  1387.       DF(J-1) = 0.0
  1388.       PNUM = FN(J-1) * XBAR(J-1) + FN(J) * XBAR(J)
  1389.       PDEN = FN(J-1) + FN(J)
  1390.       PX=PNUM/PDEN
  1391. C
  1392. C              IS IT NECESSARY TO GO FURTHER BACK?
  1393. C
  1394.       JS=J-1
  1395.       IF (JS.EQ.1) GO TO 30
  1396.       JM2=J-2
  1397.       DO 29 I=1,JM2
  1398.       IF (XBAR(JS-1).LE.PX) GO TO 30
  1399.       JS=JS-1
  1400.       DF(JS) = 0.0
  1401.       PNUM=PNUM+XBAR(JS)*FN(JS)
  1402.       PDEN=PDEN+FN(JS)
  1403.       PX=PNUM/PDEN
  1404.    29 CONTINUE
  1405.    30 DO 40 J1=JS,J
  1406.       XBAR(J1)=PX
  1407.    40 CONTINUE
  1408.    50 CONTINUE
  1409. C
  1410. C              CALCULATE STATISTIC
  1411. C
  1412.       ODF = -1.0
  1413.       OBSS = 0.0
  1414.       DO 60 I = 1, NSAM
  1415.       ODF = ODF + DF(I)
  1416.       OBSS = OBSS + FN(I) * (XBAR(I) - TBAR)**2
  1417.    60 CONTINUE
  1418.       OWSS = BSS + WSS - OBSS
  1419. C
  1420.       FRAT = (BSS / TDF) / (WSS / (TN - TDF + 1.0))
  1421.       PF = FDFC(FRAT,TDF,TN-TDF+1.0)
  1422.       OFRAT = 0.0
  1423.       IF (ODF.GT.0.0)
  1424.      *   OFRAT = (OBSS / ODF) / (OWSS / (TN - ODF + 1.0))
  1425.       CALL PBART(OFRAT,PVAL,NSAM,TN-ODF+1.0,1,SE)
  1426. C
  1427. C        PUT BOUNDS ON THE LACK OF FIT
  1428. C
  1429.       DIFFU = ((BSS - OBSS) / TDF) / (WSS / (TN - TDF + 1.0))
  1430.       PDIFFU = FDFC(DIFFU,TDF,TN-TDF+1.0)
  1431.       DIFFL = ((BSS - OBSS) / 1.0) / (WSS / (TN - TDF + 1.0))
  1432.       PDIFFL = 0.5 * FDFC(DIFFL,1.0,TN-TDF+1.0)
  1433.       WRITE (IOUT,70) BSS,OBSS,BSS-OBSS,WSS
  1434.    70 FORMAT (/' Sums of squares:',10X,'  between groups  ',G12.5
  1435.      *       //10X,'                monotone increase  ',G12.5
  1436.      *        /10X,'lack of fit  (between - monotone)  ',G12.5
  1437.      *       //10X,'                     within group  ',G12.5)
  1438. C
  1439.       WRITE (IOUT,72) FRAT, PF
  1440.    72 FORMAT(/' Overall F-ratio =',F8.2,5X,'P-value =',F6.3)
  1441.       WRITE (IOUT,80) OFRAT
  1442.    80 FORMAT(/'         F-ratio for order restriction = ',F8.2)
  1443.       IF (NSAM.EQ.3 .OR.NSAM.EQ.4) WRITE(IOUT,90) PVAL
  1444.    90 FORMAT ('                       P-value (exact) = ',F8.3)
  1445.       IF (NSAM.NE.3 .AND.NSAM.NE.4) WRITE(IOUT,100) PVAL
  1446.   100 FORMAT (' P-value (assuming equal weights) = ',F8.3)
  1447.       WRITE(IOUT,110) PDIFFL, PDIFFU
  1448.   110 FORMAT(/'     Lack of fit :  ',F6.3,' <= P-value <=',F6.3)
  1449. C
  1450.       RETURN
  1451.       END
  1452. C
  1453.       SUBROUTINE PBART(STAT,PVAL,IK,DFD,IFC,WT)
  1454. C
  1455. C        P-VALUE FOR BARTHOLOMEW'S TEST IN THE CASE OF EQUAL 
  1456. C        SAMPLE SIZES (ARBITRARY FOR IK = 3,4)
  1457. C
  1458. C        IK -- NUMBER OF GROUPS
  1459. C        DFD -- DENOMINATOR DF FOR F; DUMMY FOR CHI SQUARE
  1460. C        IFC = 1, F DISTRIBUTION
  1461. C              2, CHI SQUARE DISTRIBUTION
  1462. C        WT -- WEIGHTS:  COL TOTS FOR CHI SQUARE
  1463. C                        SE'S FOR F
  1464. C
  1465. C        CALLING PROGRAMS HAVE VERIFIED THAT THERE ARE FEWER 
  1466. C        THAN 10 GROUPS
  1467. C
  1468.       DIMENSION P(10,10), WT(IK)
  1469.       DATA PI /3.14159265/
  1470. C
  1471.       IF (IK.EQ.3) GOTO 150
  1472.       IF (IK.EQ.4) GOTO 180
  1473. C
  1474.       P(1,1) = 1.0
  1475.       DO 10 K = 2, IK
  1476.       P(1,K) = 1.0 / FLOAT(K)
  1477.       P(K,K) = P(K-1,K-1) * P(1,K)
  1478.    10 CONTINUE
  1479. C
  1480.       DO 20 K = 3, IK
  1481.       DO 20 L = 2, K-1
  1482.       P(L,K) = (P(L-1,K-1)+ FLOAT(K-1) * P(L,K-1)) / FLOAT(K)
  1483.    20 CONTINUE
  1484. C
  1485.       IF (IFC.EQ.2) GOTO 100
  1486. C
  1487. C        F DISTRIBUTION
  1488. C
  1489.       PVAL = 0.0
  1490.       DO 30 L = 2, IK
  1491.    30 PVAL = PVAL + P(L,IK) * FDFC(STAT,FLOAT(L-1),DFD)
  1492.       RETURN
  1493. C
  1494.   100 PVAL = 0.0
  1495.       DO 110 L = 2, IK
  1496.   110 PVAL = PVAL + P(L,IK) * CHIDFC(STAT,L-1)
  1497.       RETURN
  1498. C
  1499. C        IK = 3  (SPECIAL CASE)
  1500. C
  1501.   150 RHO12 = - SQRT(WT(1) * WT(3) / 
  1502.      *   ((WT(1) + WT(2)) * (WT(2) + WT(3))))
  1503.       IF (IFC.EQ.1) PVAL = 
  1504.      *   (1.0 - ACOS(RHO12) / PI) * FDFC(STAT,2.0,DFD)
  1505.      *   + FDFC(STAT,1.0,DFD)
  1506.       IF (IFC.EQ.2) PVAL = (1.0 - ACOS(RHO12) / PI) * CHIDFC(STAT,2)
  1507.      *   + CHIDFC(STAT,1)
  1508.       PVAL = 0.5 * PVAL
  1509.       RETURN
  1510. C
  1511. C        IK = 4  (SPECIAL CASE)
  1512. C
  1513.   180 RHO12 = - SQRT(WT(1) * WT(3) / 
  1514.      *   ((WT(1) + WT(2)) * (WT(2) + WT(3))))
  1515.       RHO23 = - SQRT(WT(2) * WT(4) / 
  1516.      *   ((WT(2) + WT(3)) * (WT(3) + WT(4))))
  1517.       RHO231 = - SQRT((WT(1) + WT(2)) * WT(4) / 
  1518.      *   ((WT(1) + WT(2) + WT(3)) * (WT(3) + WT(4))))
  1519.       RHO132 = - SQRT(WT(1) * WT(4) / 
  1520.      *   ((WT(1) + WT(2) + WT(3)) * (WT(2) + WT(3) + WT(4))))
  1521.       RHO123 = - SQRT(WT(1) * (WT(3) + WT(4)) / 
  1522.      *   ((WT(1) + WT(2)) * (WT(3) + WT(4))))
  1523. C
  1524.       P(2,4) = 0.125 + (ACOS(RHO12) + ACOS(RHO23))/ (4.0 * PI)
  1525.       P(3,4) = 0.75 - (ACOS(RHO231) + ACOS(RHO132) +ACOS(RHO123))
  1526.      *   / (4.0 * PI)
  1527.       P(4,4) = 0.50 - P(2,4)
  1528.       PVAL = 0.0
  1529.       DO 200 L = 2, 4
  1530.       IF (IFC.EQ.1)
  1531.      *   PVAL = PVAL + P(L,4) * FDFC(STAT,FLOAT(L-1),DFD)
  1532.       IF (IFC.EQ.2)
  1533.      *   PVAL = PVAL + P(L,4) * CHIDFC(STAT,L-1)  
  1534.   200 CONTINUE
  1535.       RETURN
  1536.       END
  1537. C
  1538.       SUBROUTINE BINOM(IIN,IOUT)
  1539. C
  1540. C        BINOMAL DISTRIBUTION
  1541. C
  1542.    10 WRITE(IOUT,20)
  1543.    20 FORMAT(/' Enter n, p, k:  '$)
  1544.       READ(IIN,*,ERR=10) XN,P,XK
  1545.       K = XK
  1546. C
  1547.       IF (P.LE.0.0 .OR. P.GE.1.0 .OR. FLOAT(K).NE.XK
  1548.      *   .OR. XK.GT.XN) GOTO 10
  1549. C
  1550.       PLO = 0.0
  1551.       XN1 = XN + 1.0
  1552.       DO 70 I = 0, K
  1553.       XI = I
  1554.       TERM = ALGAMA(XN1) - ALGAMA(XI+1.0) - 
  1555.      *   ALGAMA(XN1-XI) + XI * ALOG(P) + (XN-XI) * ALOG(1.0-P)
  1556.       IF (TERM.GT.-78.0) PLO = PLO + EXP(TERM)
  1557.    70 CONTINUE
  1558.       IF (TERM.GT.-78.0) TERM = EXP(TERM)
  1559.       IF (TERM.LT.0.0) TERM = 0.0
  1560.       PHI = 1.0 - PLO + TERM
  1561. C
  1562.       WRITE(IOUT,80) IFIX(XN),P,K,TERM,K,PLO,K,PHI
  1563.    80 FORMAT(/' Prob (binomial(',I3,',',F5.3,')  =',I3,')  = ',F7.4/
  1564.      *   27X,'<=',I3,')  = ',F7.4/
  1565.      *   27X,'>=',I3,')  = ',F7.4)
  1566. C
  1567.       RETURN
  1568.       END
  1569. C
  1570.       SUBROUTINE POISN(IIN,IOUT)
  1571. C
  1572. C        POISSON DISTRIBUTION
  1573. C
  1574.    10 WRITE(IOUT,20)
  1575.    20 FORMAT(/' Enter mean and quantile:  '$)
  1576.       READ(IIN,*,ERR=10) XLAM,XK
  1577.       K = XK
  1578. C
  1579.       IF (XLAM.LT.0.0 .OR. FLOAT(K).NE.XK) GOTO 10
  1580. C
  1581.       PLO = 0.0
  1582.       DO 70 I = 0, K
  1583.       XI = I
  1584.       TERM = -XLAM + XI * ALOG(XLAM) - ALGAMA(XI+1.0)
  1585.       IF (TERM.GT.-78.0) PLO = PLO + EXP(TERM)
  1586.    70 CONTINUE
  1587.       IF (TERM.GT.-78.0) TERM = EXP(TERM)
  1588.       IF (TERM.LT.0.0) TERM = 0.0
  1589.       PHI = 1.0 - PLO + TERM
  1590. C
  1591.       WRITE(IOUT,80) XLAM,K,TERM,K,PLO,K,PHI
  1592.    80 FORMAT(/' Prob (Poisson(',G10.4,')  =',I3,')  = ',F7.4/
  1593.      *   27X,'<=',I3,')  = ',F7.4/
  1594.      *   27X,'>=',I3,')  = ',F7.4)
  1595. C
  1596.       RETURN
  1597.       END
  1598. C
  1599.       SUBROUTINE FISHER(IIN,IOUT)
  1600. C
  1601. C        FISHER'S EXACT TEST FOR 2*2 CONTINGENCY TABLES
  1602. C
  1603.       HYPR(X11,RT1,RT2,CT1,XN) = ALGAMA(RT1+1.0) - ALGAMA(X11+1.0)
  1604.      *   - ALGAMA(RT1-X11+1.0) + ALGAMA(RT2+1.0) - ALGAMA(CT1-X11+1.0)
  1605.      *   - ALGAMA(RT2-CT1+X11+1.0)
  1606.   110 WRITE (IOUT,120)
  1607.   120 FORMAT(/' Enter row 1:  ',$)
  1608.       READ(IIN,*,ERR=110) TABL11, TABL12
  1609.   130 WRITE (IOUT,140)
  1610.   140 FORMAT(' Enter row 2:  ',$)
  1611.       READ(IIN,*,ERR=130) TABL21, TABL22
  1612. C
  1613. C        CALCULATE ROW AND COLUMN TOTALS
  1614. C        PERMUTE TABLE SO THAT FIRST ROW AND COLUMN 
  1615. C        HAVE SMALLER TOTALS
  1616. C
  1617.       RT1 = TABL11 + TABL12
  1618.       RT2 = TABL21 + TABL22
  1619.       IF (RT1.LE.RT2) GOTO 200
  1620.       DUM = RT1
  1621.       RT1 = RT2
  1622.       RT2 = DUM
  1623.       DUM = TABL11
  1624.       TABL11 = TABL21
  1625.       TABL21  = DUM
  1626.       DUM = TABL12
  1627.       TABL12 = TABL22
  1628.       TABL22 = DUM
  1629. C   
  1630.   200 CT1 = TABL11 + TABL21
  1631.       CT2 = TABL12 + TABL22
  1632.       IF (CT1.LE.CT2) GOTO 210
  1633.       DUM = CT1
  1634.       CT1 = CT2
  1635.       CT2 = DUM
  1636.       DUM = TABL11
  1637.       TABL11 = TABL12
  1638.       TABL12  = DUM
  1639.       DUM = TABL21
  1640.       TABL21 = TABL22
  1641.       TABL22 = DUM
  1642. C   
  1643.   210 GRAND = RT1 + RT2
  1644.       HYPR0 = ALGAMA(CT1+1.0) + ALGAMA(GRAND-CT1+1.0)
  1645.      *   - ALGAMA(GRAND+1.0)
  1646.       PVALL = 0.0
  1647.       PVALU = 0.0
  1648.       TERM = HYPR0 + HYPR(TABL11,RT1,RT2,CT1,GRAND) 
  1649.       IF (TERM.LT.-78.0) GO TO 280
  1650.       TERMO = EXP(TERM)
  1651.       IMAX = MIN1(RT1,CT1)
  1652.       I11 = TABL11
  1653. C
  1654.       PVALL = 0.0
  1655.       IENDL = -1
  1656.       DO 250 I = 0, I11
  1657.       TERM = HYPR0 + HYPR(FLOAT(I),RT1,RT2,CT1,GRAND) 
  1658.       IF (TERM.LT.-78.0) GOTO 250
  1659.       TERM = EXP(TERM)
  1660.       IF (TERM.GT.TERMO) GOTO 260
  1661.       IENDL = I
  1662.       PVALL = PVALL + TERM
  1663.   250 CONTINUE
  1664.       PVAL1 = PVALL
  1665. C
  1666.   260 PVALU = 0.0
  1667.       IENDU = -2
  1668.       DO 270 I = IMAX, I11, -1
  1669.       TERM = HYPR0 + HYPR(FLOAT(I),RT1,RT2,CT1,GRAND) 
  1670.       IF (TERM.LT.-78.0) GOTO 270
  1671.       TERM = EXP(TERM)
  1672.       IF (TERM.GT.TERMO) GOTO 280
  1673.       IENDU = I
  1674.       PVALU = PVALU + TERM
  1675.   270 CONTINUE
  1676.       PVAL1 = PVALU
  1677. C
  1678.   280 PVAL2 = PVALL + PVALU
  1679.       IF (IENDL.NE.IENDU) GOTO 290
  1680.       PVAL1 = AMIN1(PVALL, PVALU)
  1681.       PVAL2 = 1.0
  1682.   290 WRITE(IOUT,300) PVAL1, PVAL2
  1683.   300 FORMAT(/' Fisher''s exact test:   1 tailed  ',F7.4/
  1684.      *         '                        2 tailed  ',F7.4)
  1685. C
  1686.       RETURN
  1687.       END
  1688. C
  1689.       FUNCTION XINBTA(P, Q, ALPHA)
  1690. C
  1691. C        ALGORITHM AS 109 APPL.STATIST. (1977), VOL.26, NO.1
  1692. C        (REPLACING ALGORITHM AS 64  APPL. STATIST. (1973),
  1693. C        VOL.22, NO.3)
  1694. C        ***[FAULT INDICATOR REMOVED BY G.E.D.]***
  1695. C
  1696. C        COMPUTES INVERSE OF INCOMPLETE BETA FUNCTION
  1697. C        RATIO FOR GIVEN POSITIVE VALUES OF THE ARGUMENTS
  1698. C        P AND Q, ALPHA BETWEEN ZERO AND ONE.
  1699. C        LOG OF COMPLETE BETA FUNCTION, BETA, IS ASSUMED TO BE KNOWN.
  1700. C        ***[CALCULATION OF BETA ADDED BY G.E.D]***
  1701. C
  1702.       LOGICAL INDEX
  1703. C
  1704. C        DEFINE ACCURACY AND INITIALIZE
  1705. C
  1706.       DATA ACU /1.0E-14/
  1707.       XINBTA = ALPHA
  1708. C
  1709. C        TEST FOR ADMISIBILITY OF PARAMETERS
  1710. C
  1711.       IF (P.LE.0.0.OR.Q.LE.0.0) STOP 50
  1712.       IF (ALPHA.LT.0.0.OR.ALPHA.GT.1.0) STOP 51
  1713.       IF (ALPHA.EQ.0.0.OR.ALPHA.EQ.1.0) RETURN
  1714. C
  1715.       BETA=ALGAMA(P)+ALGAMA(Q)-ALGAMA(P+Q)
  1716. C
  1717. C        CHANGE TAIL IF NECESSARY
  1718. C
  1719.       IF (ALPHA .LE. 0.5) GOTO 1
  1720.       A = 1.0 - ALPHA
  1721.       PP = Q
  1722.       QQ = P
  1723.       INDEX = .TRUE.
  1724.       GOTO 2
  1725.     1 A = ALPHA
  1726.       PP = P 
  1727.       QQ = Q
  1728.       INDEX = .FALSE.
  1729. C
  1730. C        CALCULATE INITIAL APPROXIMATION
  1731. C
  1732.     2 R = SQRT(-ALOG(A * A))
  1733.       Y = R - (2.30753 + 0.27061 * R) / 
  1734.      *             (1.0 + (0.99229 + 0.04481 * R) * R)
  1735.       IF (PP .GT. 1.0 .AND. QQ .GT. 1.0) GOTO 5
  1736.       R = QQ +QQ
  1737.       T = 1.0 / (9.0 * QQ)
  1738.       T = R * (1.0 - T + Y * SQRT(T)) ** 3
  1739.       IF (T .LE. 0.0) GOTO 3
  1740.       T = (4.0 * PP + R - 2.0) / T
  1741.       IF (T .LE. 1.0) GOTO 4
  1742.       XINBTA = 1.0 - 2.0 / (T + 1.0)
  1743.       GOTO 6
  1744.     3 XINBTA = 1.0 - EXP((ALOG((1.0 - A) * QQ) + BETA) / QQ)
  1745.       GOTO 6
  1746.     4 XINBTA = EXP((ALOG(A * PP) + BETA) / PP)
  1747.       GOTO 6
  1748.     5 R = (Y * Y - 3.0) / 6.0
  1749.       S = 1.0 / (PP + PP - 1.0)
  1750.       T = 1.0 / (QQ + QQ - 1.0)
  1751.       H = 2.0 / (S + T)
  1752.       W = Y * SQRT(H + R) / H - (T - S) * 
  1753.      *   (R + 5.0 / 6.0 - 2.0 / (3.0 * H))
  1754.       XINBTA = PP / (PP + QQ * EXP(W + W))
  1755. C
  1756. C        SOLVE FOR X BY A MODIFIED NEWTON-RAPHSON METHOD,
  1757. C        USING THE FUNCTION BETAIN.
  1758. C
  1759.     6 R = 1.0 - PP
  1760.       T = 1.0 - QQ
  1761.       YPREV = 0.0
  1762.       SQ = 1.0
  1763.       PREV = 1.0
  1764.       IF (XINBTA .LT. 0.0001) XINBTA = 0.0001
  1765.       IF (XINBTA .GT. 0.9999) XINBTA = 0.9999
  1766.     7 Y = BETAIN(XINBTA, PP, QQ)
  1767.       Y = (Y - A) * EXP(BETA +
  1768.      *             R * ALOG(XINBTA) + T * ALOG(1.0 - XINBTA))
  1769.       IF (Y * YPREV .LE. 0.0) PREV = SQ
  1770.       G = 1.0
  1771.     9 ADJ = G * Y
  1772.       SQ = ADJ * ADJ
  1773.       IF (SQ .GE. PREV) GOTO 10
  1774.       TX = XINBTA - ADJ
  1775.       IF (TX .GE. 0.0 .AND. TX .LE. 1.0) GOTO 11
  1776.    10 G = G / 3.0
  1777.       GOTO 9
  1778.    11 IF (PREV .LE. ACU) GOTO 12
  1779.       IF (Y * Y .LE. ACU) GOTO 12
  1780.       IF (TX .EQ. 0.0 .OR. TX .EQ. 1.0) GOTO 10
  1781.       IF (TX .EQ. XINBTA) GOTO 12
  1782.       XINBTA = TX
  1783.       YPREV = Y
  1784.       GOTO 7
  1785.    12 IF (INDEX) XINBTA = 1.0 - XINBTA
  1786.       RETURN
  1787.       END
  1788. C
  1789.       SUBROUTINE MANTEL(IIN,IOUT)
  1790. C
  1791. C        MANTEL-HAENSZEL ESTIMATE OF COMMON ODDS RATIO
  1792. C        MANTEL-HAENSZEL STATISTIC TESTING SIGNIFICANCE OF 
  1793. C          COMMON ODDS RATIO
  1794. C        APPROXIMATE C.I. FOR COMMON ODDS RATIO BASED ON 
  1795. C          FLEISS(1981, EQ. 10.24).  THIS FORMULA CAN ALSO
  1796. C          BE USED TO CONSTRUCT AN APPROXIMATE C.I. FOR A 
  1797. C          SINGLE 2*2 TABLE.
  1798. C
  1799.       LOGICAL LZERO
  1800.       LZERO = .FALSE.
  1801.    90 WRITE(IOUT,100)
  1802.   100 FORMAT(/' Enter the number of 2 * 2 tables:  '$)
  1803.       READ(IIN,*,ERR=90) NTAB
  1804.       IF (NTAB.LT.1 .OR. NTAB.GT.6) GOTO 90
  1805. C
  1806.       SN = 0.0
  1807.       SD = 0.0
  1808.       SNMH = 0.0
  1809.       SDMH = 0.0
  1810.       ODDSN = 0.0
  1811.       ODDSD = 0.0
  1812.       HOMO = 0.0
  1813. C
  1814.       DO 300 K = 1, NTAB
  1815.   110 WRITE (IOUT,120) K
  1816.   120 FORMAT(/' Enter row 1 of table',I2,':  '$)
  1817.       READ(IIN,*,ERR=110) TAB11, TAB12
  1818.   130 WRITE (IOUT,140) K
  1819.   140 FORMAT(' Enter row 2 of table',I2,':  '$)
  1820.       READ(IIN,*,ERR=130) TAB21, TAB22
  1821.       IF (TAB11 * TAB12 * TAB21 * TAB22 .EQ. 0.0) LZERO = .TRUE.
  1822.       SUM = TAB11 + TAB12 + TAB21 + TAB22
  1823. C
  1824. C        MANTEL-HAENSZEL CHI-SQUARE STATISTIC
  1825. C
  1826.       SN = SN + TAB11 - (TAB11 + TAB12) *
  1827.      *   (TAB11 + TAB21) / SUM
  1828.       SD = SD + (TAB11 + TAB12) * (TAB11 + TAB21)
  1829.      *   * (TAB21 + TAB22) * (TAB12 + TAB22) /
  1830.      *   (SUM * SUM * (SUM-1.0))
  1831. C
  1832. C        MANTEL-HAENSZEL ESTIMATE OF COMMON ODDS RATIO
  1833. C
  1834.       SNMH = SNMH + TAB11 * TAB22 / SUM
  1835.       SDMH = SDMH + TAB12 * TAB21 / SUM
  1836. C
  1837. C        C.I. FOR ODDS RATIO BASED ON COMBINING LOG ODDS (FLIESS, 10.20)
  1838. C
  1839.       IF (LZERO) GOTO 300
  1840.       WEIGHT = 1.0 / TAB11 + 1.0 / TAB12 + 1.0 / TAB21 + 1.0 / TAB22
  1841.       WEIGHT = 1.0 / WEIGHT
  1842.       ODDSD = ODDSD + WEIGHT
  1843.       ORATL = ALOG(TAB11 * TAB22 / (TAB12 * TAB21))
  1844.       ODDSN = ODDSN + WEIGHT * ORATL
  1845.       HOMO = HOMO + WEIGHT * ORATL**2
  1846. C
  1847.   300 CONTINUE
  1848. C
  1849. C        MANTEL-HAENSZEL CHI-SQUARE
  1850. C
  1851.       CHISQC = (ABS(SN) - 0.5)**2 / SD
  1852.       PVALC = 2.0 * ALNORM(SQRT(CHISQC),.TRUE.)
  1853.       CHISQ = SN**2 / SD
  1854.       PVAL = 2.0 * ALNORM(SQRT(CHISQ),.TRUE.)
  1855. C
  1856. C        MANTEL-HAENSZEL ESTIMATE
  1857. C
  1858.       IF (SNMH * SDMH .NE. 0.0) GOTO 350
  1859.       WRITE(IOUT,340)
  1860.   340 FORMAT(/' Can''t evaluate odds ratio with counts of 0',
  1861.      *   ' in the summary table.')
  1862.       RETURN
  1863. C
  1864.   350 ODDS = SNMH / SDMH
  1865. C
  1866.       WRITE(IOUT,360) ODDS,CHISQ,PVAL,CHISQC,PVALC
  1867.   360 FORMAT(/' Common odds ratio = ',G12.5//
  1868.      *   ' Mantel-Haenszel chi-square statistic (P-value) '/
  1869.      *   '           uncorrected:  ',G11.4,'  (',F6.4,')'/
  1870.      *   '  continuity corrected:  ',G11.4,'  (',F6.4,')'/)
  1871. C
  1872. C        C.I. FOR COMMON ODDS RATIO
  1873. C
  1874.       IF (LZERO) RETURN
  1875.   410 WRITE(IOUT,420)
  1876.   420 FORMAT(/' Enter desired level of confidence (0,1):  '$)
  1877.       READ(IIN,*,ERR=410) CONF
  1878.       COEF = -GAUINV((1.0 - CONF) / 2.0, IFAULT)
  1879.       CONF = 100.0 * CONF
  1880.       ODDSC = ODDSN / ODDSD
  1881.       ODDSSE = 1.0 / SQRT(ODDSD)
  1882.       CL = EXP(ODDSC - COEF * ODDSSE)
  1883.       CU = EXP(ODDSC + COEF * ODDSSE)
  1884. C
  1885.       WRITE(IOUT,430) CONF,CL,CU
  1886.   430 FORMAT(
  1887.      *   ' Approximate ',F4.1,'-% confidence interval for the odds',
  1888.      *   ' ratio:'/'      (',F7.3,',',F7.3,')')
  1889. C
  1890. C        HOMOGENEITY OF ODDS RATIOS 
  1891. C
  1892.       IF (NTAB.EQ.1) RETURN
  1893.       HOMO = HOMO - ODDSN**2 / ODDSD
  1894.       PHOMO  = CHIDFC(HOMO,NTAB-1)
  1895.       WRITE(IOUT,440) HOMO,PHOMO
  1896.   440 FORMAT(/' Homogeneity of the odds ratios:'/
  1897.      *   '     Chi-square statistic = ',G11.4/
  1898.      *   '     P-value = ',F6.4/)
  1899. C
  1900.       RETURN
  1901.       END
  1902.